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A DIRECT-INVERSE METHOD FOR TRANSONIC AND SEPARATED 
FLOWS ABOUT AIRFOILS 

Leland A. Carlson 
Aerospace Engineering Department 
Texas A&H University 

SUMMARY 

A direct-inverse technique and computer program called TAMSSP that can be 
used for the analysis of the flow about airfoils at subsonic and low transonic 
freestream velocities is presented. The method is based upon a direct-inverse 
nonconservative full potential inviscid method* a Thwaites laminar boundary layer 
technique* and the Barnwell turbulent momentum integral scheme; and it is 
formulated using Cartesian coordinates. Since the method utilizes inverse 
boundary conditions in regions of separated flow* it is suitable for predicting the 
flowfield about airfoils having trailing edge separated flow under high lift 
conditions. Comparisons with experimental data indicate that the method should 
be a useful tool for applied aerodynamic analyses. 
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SYMBOLS 


isentropic speed of sound 

boundary layer coefficient in separated pressure correlation 
pressure coefficient 
Mach number 
velocity 

velocity component in the x or y direction respectively 
transformed velocity at boundary layer edge 
velocity in the boundary layer 

law-of-the-wall and law-of-the-wake velocity parameters 
Cartesian coordinates 
angle of attack 

ratio of specific heats* assumed to be 1.4 
circulation 

boundary layer thickrass 
polar coordinate 
computational coordinates 
potential function 
perturbation potential 



SYMBOLS (Continued) 


Subscripts: 

c* fretstroam condition 

b body 

• boundary layer edge 

i»J grid location 

LB folding edge 

TS trailing edge 

J't'Xiy differentiation 
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GENERAL DESCRIPTION 


INTRODUCTION 

Over the past decide * sever*! finite-difference potential flow methods 1-3 
have been developed and successfully used for the design and analysis of 
subsonic and transonic airfoils at and near cruise conditions. However in the 
analysis of high performance airfoils* aerodynamicists would also like to be able 
to predict airfoil pressure distributions and aerodynamic coefficients at high lift* 
high angle-of-attack conditions. Since such situations are frequently 
characterised by regions of separated flow on the upper surface and are 
dominated by strong viscous interaction effects* inviscid methods alone are not 
applicable. Furthermore* subsonic-transonic analysis method»3*4 W hich couple 
inviscid and boundary layer solutions typically only include the effects of weak 
viscous interaction and generally fail to give accurate results when separated 
flow exists on the upper surface. 

However* it has been demonstrated® - ® that the direct-inverse technique 
coupled to a suitable boundary layer method can be successfully applied to low 
speed flows about airfoils having massive separation. In addition* Barnwell?* 
Dvorak and Choi 10 , and Tavern* 11 have developed similar methods for transonic 
flows. Barnwell's method* however* is limited in application in that it utilizes 
for its inviscid solver the transonic small perturbation equation. Further* 
References 9 and il only include the effects of viscous interaction due to a 
turbulent boundary layer. On most airfoils* particularly on the lower surface at 
high angles of attack, extensive regions of laminar flow exist. 
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This report describee a flow model and computer program* called TAHSEP, 
which can be used to predict the flowfield about a single element transonic air-foil 
at high angle of attack high lift conditions with trailing edge separation. Since 
the method is based upon the TRANDES 4 ind TRANSEP* codes, it can also be 
used for subsonic-transonic analyses not Involving separation. 

METHOD OF APPROACH 

The present approach is based upon the direct-inverse method developed in 
the TRANDBS and TRANSSP programs and the ability of this method to use either 
the displacement surface (airfoil ordinate plus displacement thickness) or 
pressure as the airfoil boundary condition. For the high angle-of-attack case, 
the airfoil lower surface only experiences weak viscous interaction and 
frequently has a long laminar run before transitioning to fully turbulent flow. 
Thus* the present model includes an initial laminar boundary layer calculation in 
its viscous mterction section. On the upper surface the boundary layer is also 
intially laminar* but it quickly becomes turbulent in character followed m many 
cases by boundary layer separation and a separated zone which can extend over a 
significant portion of the airfoil surface. In the present model* this separated 
region is treated inversely in that the pressure distribution along the effective 
displacment surface streamline is determined iteratively as part of the solution 
and used as the airfoil boundary condition. Consequently* the present method has 
been modeled as shown on Figure i. 

To obtain the inviscid portion of the flowtield* the full potential equation for 
two-dimensional compressible flow 13 used in nonconservattve form as 

-aUJ.y *(•'• },•)!„ ’O 
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wham the subscripts dsnott partial differentiation. By defining a perturbation 
potential* f, such that 

& s T Y <)*>*"'* +*)<*$ (2) 

where the velocity components are given by 

U* 7-< coM '<&) 

(3) 

V - §y : +tP r ) 

the governing equation in terms of the perturbation potential can be written as 
(a. 1 -O') <P„ -ZU\/d f tU‘-l/‘)<Py r *0 HI 


with 

a 4j: ^ ~1+ ) «> 


The nonconservative form of the potential equation was selected for the 
present problem because for two-dimensional flows results obtained with it 
agree better with Eule^ solutions than those obtained using the fully 
conservative form of the equation.i2 j n idditioni the conservative formulation 
appears to break down in two-dimensional cases shortly after the onset of 
supercritical flow.13 
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In the present model, Equation* (3-5) are finite differenced using a roated 
difference scheme and solved iteratively using column relaxation in a stretched 
Cartesian grid which maps the infinite domain to a finite computational box. The 
appropriate boundary condition at infinity is 

<P= To n V (/-*£( {6) 

where 9 is the polar angle* and P is the circulation* which is determined by the 
change in potential across the Kutta-Joukowski cut at the trailing edge of the 
airfoil. 


Likewise* the appropriate airfoil boundary condition in the direct regions 
(regions without separation having only weak viscous interaction) is the flow 
tangency condition given by the ordinates of the airfoil displacement surface* i.e. 


( 4A - (X) a h 

\ " \ V'k ^ ( cia * (j) xb ) 


(?) 


In the inverse or separated flow region the pressure distribution along the 
effective displacement surface streamline is considered specified and used as the 
boundary condition. As shown in Reference 2* this approach leads to a derivative 
boundary condition for the inverse region of the form 


^ 


( 8 ) 


Complete details concerning the finite difference scheme, the stretched Cartesian 
grid system, and the treatment of the boundary conditions are given m 
References 1.2* and 4. 
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To include viscous effects* the basic Approach is to calculate a boundary layer 
displacement thickness for the weak interaction regions and to use it to correct 
the location of the displacement surface (i.e., airfoil ordinate plus displacement 
thickness). For the strongly interacting separated zone on the upper surface, the 
pressure is determined from the interaction solution and the location of the 
displacement surface is computed by integrating the surface tangency condition. 
Eq. (7). with the initial conditions specified by the displacement surface 
ordinates at the separation point, which is the interface between the two regions. 
The location and slopes of the displacement surfaces are updated regularly 
throughout the iterative solution. 

In the present method, the laminar portion of the boundary layer is computed 
using a compressible Thwaites method similar to that used previously in 
TRANSEP.6 jhe transition location is determined from a Granville tjpe 
correlation^ based upon the difference between the local momentum thickness 
Reynolds number and the value at the laminar instability point combined with the 
pressure gradient history. Sometimes, particularly on the upper surface at high 
angles of attack, laminar separation is predicted upstream of the transition 
point. In these cases, the local momentum thickness Reynolds number is compared 
to an empirical correlation in order to determine if the laminar bubble is long or 
short. If the bubble is short, its length is assumed to be one horizontal delta-x 
grid width and the turbulent flow computation is initiated at the next downstream 
grid point. If the estimate indicates that the bubble is long, the calculation 
proceeds, but a warning is printed which indicates that the results are probably 
in error. 
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After transition* th* turbultnt boundary layer is computed using the 
Simplified Kuhn find Jfiolsen method (SKAN) as developed by Barnwell in Reference 
9. This method was selected because it is efficient* reliable* and yields excellent 
predictions of displacement thicknesses and separation point location. The SKAN 
turbulent boundary layer method solves the integral forms of the momentum 
equation* moment of momentum equation* and the derivative of the Coles' 
law-of-the-wall law-of-the wake relationship applied at the boundary layer edge. 
After considerable effort9 ( these equations can be transformed into a set of 
simultaneous ordinary differential equations* i.e. 



which can be solved for the wall friction velocity* u** the wake parameter u p * 
and the boundary layer thickness <f * using a second-order predictor-corrector 
technique . The remaining quantities of interest such as displacment thickness 
and momentum thickness are then determined from these variables. The numerical 
integration is terminated at the separation point* where the wall friction 
velocity* u # , vanishes. 

The method uses a two-layer eddy-viscosity model* which ignores the viscous 
sublayer terms* consisting of an inner layer Prandtl mixing length model and an 
outer layer Clauser model; and the intermittency factor as well as several 
density ratios appearing in the fundamental equations are approximated. In 
addition* the method assumes an adiabatic wall. 
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Thus, on the low*r surface the flow is computed using direct boundary 

conditions (airfoil specified) including the effects of weak viscous interction. On 

the upper surface ( the flowfield is also computed directly with viscous 

interaction up to the separation point, which is determined as part of the 

boundary layer solution. Downstream of separation, inverse boundary conditions 

are utilized, and the pressure must be specified. Fortunately, if the skin friction 

at the wall is assumed to be zero in the separated zone, the SKAN formulation can 

be used to obtain a closed form solution for the velocity, and hence the pressure, 

at the outer edge of the separated zone. The resultant analytic expressions for 

the velocity and pressure are 

jja _ JT C, ( cos* *«r) f 



As can be seen, the separated pressure depends upon the flowfield solution 
via the mviscid perturbation potentials at the separation point and the trailing 
edge, the size of the separated zone, and thru eg* the boundary layer solution at 
the separation point. In addition, this closed form solution predicts a variable 
pressure distribution for the separated region. At low freestream Mach numbers 
this variation is extremely small and is essentially constant. However, at 
freestream Mach numbers of 0.3 and above, the variation becomes signficant and 
influences the resultant flowfielr solution. This trend and separated pressure 
vacation is in accord with experimental observations and is a significant 
improvement over previous methods which assumed constant pressure in the 
separated zone regardless of flow conditions. However, since at low speeds the 
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separated pressure is essentially constant and the complexity introduced by 
equations <10—ll» nay not be warranted* the present method contains the option 
of either using a constant pressure in the separated region or the variable 
distribution determined by the closed term solution given above. 

In principle* the separated region and the wake should be accurately modeled 
with respect to physical phenomena and internal details* and this approach has 
been taken by other investigatorsi0-!5-16 a j n p rf8Bn t model* however* the 
wake region contains very tew computational prints since the coordinate system 
rapidly stretches to intinity. Thus* the wake is assumed to be inviscid with a 
constant pressure trailing edge tormed by the upper and lower displacement 
surfaces. Fortunately* extensive numerical experiments with the present and 
previous 7 models indicate that the pressure distribution and aerodynamic 
coefficients are primarily deoe^dent upon obtaining accurate predictions for the 
location of the separation point and the magnitude and variation of the separated 
pressure. Apparently* the details of the wake region are of secondary 
importance. Since the present method obtains the separation point location 
directly from the solution for the wall friction velocity* u®, and the pressure 
variation from a solution which couples the inviscid and viscous parts* it should 
yield reasonable engineering results. 

INTERACTION AND ITERATION PROCEDURE 

The iteration and interaction procedure used in TAMSEP is similar to that 
used in the low-speed program TRANSEP^, and it is outlined in schematic form on 
Figure 2. The program first reads all necessary input and initializes the 
perturbation potential at all grid points to zero. Next it computes 
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transformation factors tnd coordinate* associated with the stretched Cartesian 
grid tor the initial grid specified by the input data. Included in this process is 
the computation of all airfoil ordinates and slopes required on the computational 
grid. 

Since the initial grid is normally very coarse with a default size of 13 m 7, only 
fifty inviscid iterative cycles are computed on this grid. The calculation 
procedure used for the inviscid potential equation is the same iterative 
successive column over-relaxation scheme used in TRANDBSM, This limited 
number of cycles serves to rapidly create an approximate starting solution for 
succeeding grids. After these are completed* the grid is then halved and the 
solution interpolated onto a new grid, v.*uch has a default s:~s of 25 :: I?. 

After obtaining all necessary coordinates* stretching factors* airfoil 
ordinates and slopes* etc.* for this second grid, the method then performs fifty 
inviscid iterative cycles before considering any type of viscous interaction. 
Experience has shown that it is important to perform a limited number of inviscid 
cycles at the beginning of each new grid in order to eliminate any "problems" 
introduced by grid halving and interpolation. 

After these initial iterations* the program then checks to see whether or not 
the user desires viscous interaction to be included by examining the value of the 
variable ITACT. If viscous interaction is desired* which is specified by the 
ITACT default value of one* the program then checks to see if an initial laminar 
boundary layer is to be included (1LAM*1) or if the viscous calculations are to be 
for a turbulent boundary layer with user specified transition points (ILAM»0). 
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Upon completing the boundary layer computations tor the current flowfield 
solution« the program then calculates the ordinates and slopes of the upper and 
lower displacement surfaces. Since it only involves weaK viscous interaction! the 
lower surface computations are from the leading edge or the lower surface 
stagnation paint, whichever is further aft. to the trailing edge. However on the 
upper surface they are only from the leading t>dge to the separation point or to 
the trailing edge, whichever is less. This process involves smoothing of the 
displacement thiel«r>tss values! properly adding them to the airfoil ordinates, and 
spline fitting the resulting points. 

At this point the procedure depends upon whether or not separation has been 
detected on the upper surface. If separation does not exist prior to the last grid 
point on the airfoil upper surface, additional inviscid cycles are performed 
before returning to the viscous interaction loop. However, if separation is 
predicted, then the method must determine the pressure distribution and the 
location of the displacement surface in the separated zone. 


Exactly how the separated pressure distribution is determined depends upon 
the user specified variable KSEP. If KSEP is zero, the pressure is assumed to be 
constant m the separated zone and is computed by 


C P = 


- * ( f’v ~ O 


(12) 


While this expression is a small perturbation approximation for Cp (Sepi u6a g e 

has been found to be accurate and adequate for low speed incompressible flows. 
At freestream Hach numbers of r. 3 and higher, however, the variable separated 
pressure option, specified by a KSEP value of one, should be used. In this case, 
the pressure distribution along the displacement streamline in the separated 
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region is determined by Equations <10) end <1 i> above. Note that both approaches 
determine the separated zone pressure, which depends upon the currant solution, 
by conditions at both the separation point and at the trailing edge and not Just on 
conditions in the vicinity of separation. This result is in agreement with the 
conclusion of 0rossi7 that conditions at the downstream end of the separation 
zone partially influence the separation pressure level. 

After determining the separated pressure distribution corresponding to the 
predicted separation point and the current potential flow solution, the 
corresponding upper surface displacement surface must also be computed for the 
separated zone. When separation exists, the previous method of adding the 
computed displacement thickness to the original airfoil ordinates is inappropriate 
since the values for displacement thickness predicted by the SKAN method are 
probably inaccurate in separated regions. Instead, the present approach is to 
solve, using the current potential flow solution, the differential equation 

AM _ (yL) , s l21Ui£u 

\TrJ h ~ {0 'i ~ c««u(; <p u <J3) 

for the y-ordinates of the separated displacement surface as a function of x. 
Based upon previous studie»2,6, Equation <131 is solved using the Rungo-Kutta 
method of order four and the displacement surface ordinate at the separation 
point as the initial condition. Jn addition, m the process of solving this equation 
and tp^must be evaluated by finite differences. While several formulations 
are possible^, numerical studies indicate that accurate displacement surfaces are 


obtained using the following 



w ®Bsy^raaHrs»£ Tr . :*_ 
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(14b) 


In Equations (14), the point (i,rl> is the first ghost point below the displacement 
surface. Its value is determined as part of the inverse prssure boundary 
condition. 

Since the present process is iterative and the potential solution uses the 
separated pressure distribution as an inverse boundary condition, the solution of 
Eqs. (13) and (14) should yield upon convergence a separated zone displacement 
surface or free streamline that is compatible with the pressure distribution and 
potential flow solution. 


At this point m the iteration-interaction procedure a check is made to see if 
the solution has converged or it the maximum number of iterations tor a given 
grid size has been exceeded. It neither situation is true, ten more inviscid cycles 
with the new displacement surfaces and separated pressure distribution are 

i 

i performed prior to repeating the viscous interaction loop. If, however, either 

j condition is satisfied and the finest grid specified by the user has not been used, 

I 

I the grid is refined and the entire process shown on Figure 2 is repeated. If the 

last grid solution has been obtained, then a final output is printed and the 

, solution is finished. 

-i 
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It should be noted that th» calculations on a given grid art stopped and 
assumed to be converged when the maximum perturbation potential change is less 
than some user specified value. However* when separation is present it is usual 
for the calculations on each grid to be terminated due to the number of iterative 
cycles exceeding a maximum user specified value (particularly on computers which 
only retain seven significant digits). In those cases* the existence or degree of 
convergence can be determined by examining the variation in the number of 
supersonic points* the location of separation* and the traili'ig edge ordinate of 
the upper displacement surface. All these values are printed out every ten 
iterative cycles. If they stabili 2 e prior to the end of the computation on a given 
grid* then the results can be assumed to be converged. Normally* it is sufficient 
to perform 800 cycles on the coarse grid (25x13)* 400 on the medium grid (49x25)* 
and 400 on the fine grid ((97x49)* although occasionaly more may be needed. In 
determining convergence* it should be remembered that the present method is 
supposed to obtain a steady state solution. At angles of attack above maximum 
lift* the actual flowfield about an airfoil is usually unsteadyi7il8 a j n those 
cases* the present method probably will not converge and may enter some type of 
oscillatory behavior which appears to represent an unsteady flow pattern. 
However* the present method is not “time-accurate" and such results should only 
be viewed as indicative of the presence of significant unsteady phenomena. 
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USSR INSTRUCTIONS 


CODS DESCRIPTION 

The TAMSEP cod* consists of * main program and eighteen subroutines. The 
subroutines and their relationships are shown in a subroutine tree on Figure 3. 
The subroutine names and their functions are as follows: 

FOIL — Reads in initial airfoil shape and determines ordinates and slopes at 
computational points. 

VISACT — Computes turublent boundary layer when viscous interaction included. 

THWAIT — Computes the laminar boundary layer when viscous interaction 
included. 

F1T2 — Curve fit routine used by Thwaites method. 

VALUE — Initializes the flowfield to zero perturbation potential. 

SOLVE — Sets up the matrix coefficients used in the SLOR relaxation scheme. 

PRESS — Computes the pressure distribution on the airfoil. 

COORD — Sets up the coordinates in the computational and physical grids and 
computes the stretching factors. 

FLOW! — Solves flowfield in front of the airfoil. 

FLOW2 — Solves flowfield in the direct region above and below the airfoil. 

FLQW3 — Solves flowfield in the inverse region. 

WAKE — Solves the flowfield behind the airfoil. 

SHAPE — Computes the shape of the airfoil displacement surface in the 
separated zone. 

TRID — Tndiagonal equation solver. 

HALVE — Doubles the grid size and interpolates old values to obtain starting 
values on the new grid. 

PLOT Creates a printer plot of the Cp an( j displacement surface. 
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ARC — Determines the Arc length o f the Airfoil coordinates And splines the 
coord in Ate ■ versus Arc length. 

SPLINE — Computes a cubic spline thru a set of points. 

The TAHS8P code is written in FORTRAN IV progremming language And is 
designed for use on IBM# AMDAHLi CDC, DEC# And similar computers. In 
nonoverlAy mode it requires less then 320)000 bytes on An AmdAhl 470/V8. Using 
a FORTRAN H extended compiler At the optimizAtion level two* it needs About 17 
seconds for compilation and obtains a solution on a 97x49 grid in about 160 
seconds at a rate of around 15*000 points/second. Some slight modification to 
formats# etc. may be required to run the program on different computer systems 
or under a FORTRAN 77 compiler. 

INPUT DESCRIPTION 

The input to the TAM SEP code is read in eight separate blocks. The first one 
contains a user supplied title# while the second and third blocks specify all the 
floating point and integer parameters needed to run the program. These 
parameters are input via namelists and if not specified are assigned default 
values by the program. Blocks four and seven are optional and are only included 
when the parameter IREAD is one. They read a non-zero perturbation potential 
starting solution and an intial airfoil description. Blocks five and r ig! t are 
associated with input for the design option in the program and are only included 
when the parameter INV is one. Finally# block six contains the description of the 
airfoil under consideration. For an analysis computation# only blocks one# two# 
three and six would be included in the input stream. 
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DETAILED INPUT DESCRIPTION 
Input Block 1: Title 

Thi* block consists of a single line of input and is read by the main program. 

NTITLB Description of case. Up to GO alphanumeric characters. Appears on the 
printed output at the beginning of the results of each grid. 

Input Block 2: Floating Point Parameters 

This block of input is read by the main program via a namelist called F1NP. 

M Free stream Mach number (real variable). Default 0.5 

W Relaxation factor for subsonic points. Should be in the range 0<W<2.0. 

Default 1.7 

XI X location where the direct mode calculation procedure stops. In the 
analysis mode it should be set to 0.5 (i.e. the trailing edge). In the inverse 
(design) mode it is usually set to slightly less than the third point from the 
leading edge or larger. Default 0.5 

X2 End of the inverse region. For analysis cases set to a large number. In 
the inverse (design) case set to 0.5 ( i.e. the trailing edge). Default 10000.0 

ALP Angle of attack in degrees. Default 0.0 

EPS Subsonic damping factor to match difference equations at 3onic line if 
needed. Has no effect on the accuracy of the solution. Only affects stability and 
convergence rate. Normally it is not needed. Default 0.0 
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EPSS Supersonic damping factor for iterative stability. Has no effect on the 
accuracy of the converged solution, only on the stability and convergence rate. 

Should typically be about M 2 - 1 * where M max is the maximum local Hach 
number. Default 0.4 / 

X4 The positive X locations where the coordinate stretching changes. It 
should be near the airfoil trailing edge. Default 0.49 

S4 The positive psi value in the computational plane where the stretching 
changes. Default 2.0 

CONV Convergence criteria control value. Iterations stop when the maximum 
change m the perturbation potential between relaxation ./cles is less than CONB. 
Default l.E-05 

A1 Stretching constant for the Y direction. It can be used to control the Y 
and eta spacing near the horizontal axis. It is usually best to have the psi and 
eta spacing equal near the leading edge of the Airfoil. Default 0.246 

A2 First stretching constant for the X-direction. It is equivalent to (2/ff * 

(dx/di ) at . The value of A2 determines the horizontal step size near the 

leading and trailing edges, i.e. 

d^»x4*tLA2 4 j 2 TrA2 g<lt54)) 

2 2 UMAX -1) 

Bee Appendix A of Reference 4. Default 0.15 

A3 Second stretching constant for the x- direction. It determines the 
physical location of the vertical grid line adjacent to the grid side edge. Default 
3.87 
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RN Freestream Reynolds number based on chord length. Used only when 
viscous interaction (1TACT«1) included. Default 20.E+06 

XIBDLY The x -location at which upper surface transition is assumed to occur. 

i 

The turbulent boundary layer calculation starts at the next grid point. The 
relationship to percent chord is: 

XIBDLY • (Xchord - 50.01/100.0 

Used only if viscous interaction included (ITACT*i) and laminar boundary layer 
ignored (1LAM-0). Default -0.44 

C1R Circulation about airfoil. If an initial solution is input (1READ*1)| it 
must be the corresponding value of circulation (CIR»CL/2.0). Default 0.0 

t 

i 

CDCORR Correction to the wave drag coefficient. Because of the lack of a large 
number of points m the leading and trailing edge regions# the wave drag 

coefficient has an error associated with grid size> spacing) and lift coefficient. , 

I 

The magnitude of CDCORR must be determined by the user by empirical methods. i 

i 

Note that the correction should be different for each airfoil and grid combination. I 

t 

Default 0.0 


RDEL Relaxation parameter for the boundary layer displcement thickness. It is 
U3ed only whe viscous interaction is included (ITACT»i) and IMAX is less than or 
equal to 55. Default 0.25 

RDELFN Fine grid relaxation parameter for the boundary layer displacement 
thickness. It is used only when viscous interaction is included and IMAX is 
greater than 55. Default 0.125 
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SP Maximum value allowed for the Nash-Macdonald separation parameter 
when x<XSBP. Used only in the design case (lNV“i> when computing the boundary 
layer over the design surface. Default 0.004 

XSEP X location after which the Nash-Macdonald separation parameter can 
exceed SP. Used orly m the design case (IN V«i) when computing the boundary 
layer over the design surface. Default 0.44 

RCPB Not used. Ignore. 

CPD Not used. Ignore. 

XMON Not used. Ignore. 

XLSEP Should always be initialized to 0.5. Default 0.5 

XPC Location after which the lower surface displacedment thickness is 
required to continue decreasing once it has started to decrease. Upstream ' XPC 
the displacement thickness is required to be monotomcally increasing. For most 
aft-cambered airfoils it should be set to O.ii and for conventional airfoils it 
should be set to 0.5. Default 0.1 

XLBDLY The x-location at which lower surface transition is assumed to occur. 
Same relationship to chord as XIBDLY. Used only if viscous interaction included 
(I TACT” 1) and laminar boundary layer ignored (ILAM=0). Default -0.44 

RLAX Relaxation parameter for the separated pressure level in the constant 
separated pressure option <KSEP=0). Sometimes needed to enhance convergence. 
Used only when lTACT=i ( IMASS»i» and K5EP=0. Default i.O 
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RAD’JJ Leading edge radius of the airfoil nondimensionalized by the chord. 

Used only if ITACT-i and 1LAMM. Default 0.0159 

Input Block 2: Integer Parameters 

This block of input is read by the main program via a namelist called IINP. 

IHAX Number of vertical grid lines in the horizontal direction on the first grid. 
1*1 corresponds to upstream infinity and WMAX corresponds to downstream 
infinity. For each grid refinement IHAX is increased such that the new IMAX is 
two times the old value minus one. The limit on IHAX is 99. Default i3 

JMAX Number of horizontal grid lines in the vertical direction on the first grid. 
J*i corresponds to infinity below the airfoil and J* JMAX is infinity above the 
airfoil. The same formula and limit that apply to IMAX also apply to JMAX. 
Default 7 

IKASE An integer number describing the case being computed. It is limited to a 
maximum of six digits and is printed at the beginning of the pressure printer plot 
for each grid. Default 100 

INV Parameter determining the program mode. It should be zero for analysis 
cases and one for inverse design cases. Default 0 

MITER Maximum number of interations (complete relaxation cycles) allowed on 
the first grid. MITER is halved for each grid refinement. However» on the fourth 
grid# MITER is reset to 400. Default 1600 

NHALF Number of grid refinements. Default 2 
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ITACT Viscous interaction control parameter. It should be set to 2 ero for 
analysis cases without interaction and for design cases. It should be one tor 
analysis cases with interaction. Default 0 

ISKP2 Airfoil update control parameter for grid two. It should be zero if an 
airfoil shape update is desired on grid two every ten iterations. It should be one 
it an update is not desired until grid two solution is completed. Only used :n the 
invest design mode (INV«i>. Default 0 

ISKP3 Same as ISKP2 but for grid 3 

ISKP4 Same as ISKP2 but for grid 4 

ITERP Interpolation parameter for the design pressure distribution on grid 
four. If in the design mode the input pressure distribution for grid four is to be 
read as input data, INTERP should be zero. If it is desired to hnerarly 
interpolate the pressure distribution of grid three* it should be one. Default 0 

IREAD Starting solution control parameter. If IREAD is zero* the initial 
perturbation solution is assumed to be zero everywhere. If it is one* an initial 
solution is read as data. The latter would only normally occur when a user wished 
to restart a solution which had previously been saved. Default 0 

LP Relaxation cycle interval at which boundary layer details are prin^-i. For 
diagnostic purposes suggest 50 or 100. For normal information purposes* suggest 
a value of 200. Default 1000 (no printout) 
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IM ASS Massive separation parameter. It should be one if the massive 
separation option is desired in analysis cases and is active only if ITACT is one. 
In inverse design cases (INV«i) it should be zero. Default 1 

ILAH Boundary layer parameter. If zero* boundary is computed as if all 
turbulent with transition at X1BDLY and XLBDLY. If one, boundary layer is 
considered laminar-turbulent with natural transition. Default i 

IPRTi Print parameter. If one* perturbation potential values printed at the 
completion of each grid. Default 0 

IPRT2 Print parameter. If one* x and y velocities at each grid point printed at 
the completion of each grid. Default 0 

KSEP Separated pressure distribution parameter. If zero* pressure in 
separated zone is assumed to be constant. If one, pressure m separated zone is 
considered variable. Default 1 


Input Block 4: Starting Solution (Optional) 

This block of data is read by subroutine VALUE onl/ if the integer parameter 
IREAD has the value of one. 

P(I,J) Nondimensional perturbation potential at point I*J. Read by rows starting 
at J=*JMAX down to J=*l. Each row runs from 1=1 to I=IMAX and starts on a new 
line. Format 5E15.7 
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PB(1) Nondimensional perturbation potential at point 1 on the y»0 (iae. J« JB) 
grid line that is associated with the lower surface of the airfoil. Read from 1*1 
toI-lHAX. Format 3B 15.7 


Input Block 5: Direct Inverse Parameters (Inverse Design Mode only) 

This single line of input is read in subrout.ne COORD only when the Inverse 
design mode is active <1NV«1>. 

XliX2 Same definition as in Block 3. However, when the inverse design mode is 
active, these values are read prior to the solution of each grid. This block 
corresponds to the first grid; and, thus, should always, use Xl*0.5 :nd X2« 10000.0 
Format 2F 10.5 


Input Block 6 : Airfoil Description 

This block of data is read by subroutine FOIL and describes the airfoil used in 
the analysis mode or the starting airfoil for the inverse design mode. 

NI The number of coordinate pairs used to describe the upper surface of the 
airfoil. Maximum value limited lo 99. Format 15 

X1(I),YI(I> Coordinate pairs describing the upper surface of the airfoil. The 
leading edge corresponds to X1 = 0.0 and the trailing edge is XI»i.O. The vertical 
ordinate, YI, is nondimensionahzcd by chord. Read starting with l»i to I»NI. 
Format 8F10.4 
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DERIX» DBRIYi DERFX. DERFY Parameters describing tht leading and trailing 
*dgt of the airfoil. DERIX it dx/ds of the airfoil upper surface at tha laading 

adga. It is usually zaro. LERIY is dy/ds of tha airfoil uppar surface at tha 
laading adga and it is usually 1.0. DERFX is d 3 x/ds 3 of tha airfoil uppar surfaca 
at tha trailing adga. It is usually sufficiantly accurate to usa 0.0. DERFY is 
d 3 y/ds 3 of tha airfoil uppar surfaca at tha trailing adga. It is usually 
sufficiantly accurata to usa 0.0 . Format 4FI0.4 

NIB Tha number of coordinate pairs used to describe tha lower surfaca of tha 
airfoil. Maximum value limited to 99. Format 15 

XlB(l)iYlBU) Coordinate pairs describing tha lower surfaca of the airfoil. The 
leading edge corresponds to XI»0.0 and the trailing edge is Xl*1.0. The vertical 
ordinate. YIB. is nondimensionalizad by the chord. Read starting with 1*1 to 
! b NIB. Format 8F10.4 

DERIXB.DERIYB.DERFXBiDERFYB Parameters describing the leading and trailing 
edge of the airfoil. DERIXB is dx/ds of the airfoil lower surface at the leading 
edge. It is usually zero. DER1YB is dy/ds of the airfoil lower surface at the 
leading edge and it is usually -1.0. DERFXB is d 3 x/ds 3 of the airfoil lower 
surface at the trailing edge. It is usually sufficiently accurate to use 0.0. 
DERFYB is d 3 y/ds 3 of the airfoil lower surface at the trailing edge. It is usually 
sufficiently accurate to use 0.0. Format 4F10.4 
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Input Bioctt 7: Starting Airfoil Description (Optional) 

This block of data is read from subroutine FOIL and is only read if the 
integer input parameter IREAD is one. It effectively overwrites the information 
from input block six. 

YUU)*Y1 <1>*SLU(1)*SLLU) Values describing the airfoil on the starting grid. YU(!) 
and YL(I) are the upper and lower surface ordinates* nonditr.cnsior.alized by chord* 
at cnord location XU). SLU(I) and SLL(l) are the upper and lower surface slope 3 at 
chord location X(I). The X(I) values depend upon the size and spacing associated 
with the starting grid. The group of four values is read starting at the 1 value 
corresponding to the first point downstream the lead.. »j edge (1>ILE) and 
ending with the point just upstream of the trailing eoge (I»ITE>. .''"•'mat 5E15.7 

DUPOLDU)*DLWOLD(I> Values describing the boundary layer displacement 

thicknesses on the starting grid. DUPOLD(l) and DLWOLDU) are the upper and 

lower surface displacement thicknesses corresponding to the chord location X(I). j 

t 

These are read starting at the 1 value corresponding to the first point > 

downstream of the leading edge U=1LE) and ending with the point just upstream of 
the trailing edge (I»1TE). Format 5E15.7 


Input Block 8: Design Pressure Distribution (Inverse Design Mode only) 


This block of data consists of four sections which are only included in the 
inverse design mode (1NV=1). In that mode only the last three sections would 
usually be included. 
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Stction 1 — Starting solution design pressure distribution read by subroutine 
FOIL. This section would only be included if a design solution were being 
restarted (i.e. INV*1, IR5AD*i, and MHALF*1) and it would only affect the first 
grid considered. 

CPU(1) Upper surface inverse region pressure coefficient values for the design 
case starting with I«Ii* which is the first grid point after Xi and ending with 
!*ITB* the grid point just upstream of the trailing edge. Format 8E10.3 

CPL(l) Lower surface inverse region pressure coefficient values for the design 
case starting with 1*11, which is the first grid point after XI and ending with 
I*ITE* the grid point just upstream of the trailing edge. Format 8E10.3 

Section 2 — This section reads in the starting and ending points of the inverse 
region from subroutine COORD and the inverse design pressure distribution from 
subroutine FOIL for the second grid. Used only m the inverse design case 
<INV«1>. 

Xi*X2 XI is the location where the direct calculation stops and the inverse 
calculation begins. Typically* it is slightly less than the third point from the 
leading edge or larger. X2 is the location where the inverse calculation stops. It 
should always be set to 0.5. Format 2F10.5 

CPU(I) Upper surface inverse region pressure coefficient values for the design 
case starting with I*Ii* which is the first grid point after XI and ending with 
I*ITEi the grid point just upstream of the trailing edge. Format 8E10.3 
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CPL(l) Lower surf Act inverse region pressure coefficient values for the design 
case starting with I*Ii, which is the first grid point after XI and ending with 
I>ITE, the grid point Just upstream of the trailing edge. Format 8B10.3 

Section 3 — This section reads in the starting and ending points of the inverse 
region from subroutine COORD and the inverse design pressure distribution from 
subroutine FOIL for the third grid. Used only in the inverse design case (INV«i). 
The input variables and descriptions are the same as Section 2 above. 

Section 4 — This section reads in the starting and ending points of the inverse 
region from subroutine COORD and the inverse design pressure distribution from 
subroutine FOIL for the fourth grid. Used only in tne invar c design case (INV»1) 
when 1TERP is zero. Note that in References <2> and (4), the use of grid four for 
inverse design is not recommended. The input variables and descriptions are the 
same as Section 2 above. 


OUTPUT DESCRIPTION 

THe printed output when the program is operated in the inverse design mode 
is identical to that described in Reference (4). When the program is operated in 
the analysis mode with the massive separation option, the output for rrids two, 
three, and four has the form shown below. Since the first grid assumes mviscid 
flow only, its printout only includes those portions associated with an mviscid 
solution. 


30 




*? ' 

* > 

r 


i 


1. Heading (user supplied) 

2. Case Number 

3. Mach number and angle of attack 

4. Case type callouts* i.e.* 

IN VISCID ANALYSIS CASS 

WITH LAMINAR TURBULENT VISCOUS INTERACTION 
AND MASSIVE SEPARATION 

AND VARIABLE PRESSURE IN SEPARATED REGION 

5. Input data in namelists FINP and IINP 

6. Coordinate System for the current grid printed as I* X(I) followed by J, Y(J) 

7. Ordinates of the current airfoil displacement surface 

X — Horizontal ordinate* where -0.5 is the leading edge and 0.5 is the 
trailing edge 

YU — Upper displacement surface ordinate 
YL — Lower displacement surface ordinate 
UPPER SLOPE — Slope of upper displacement surface 
LOWER SLOPE — Slope of lower displacement surfcace 

8. Iteration history at ten-cycle intervals 

CIR — Circulation 

DPM — Maximum (0 correction (absolute value) in the last relaxation cycle 
with the corresponding (I,J) grid location 

NSSP — Number of supersonic grid points 

DELTAY OR YUTE — In the design case* the change in YU(ITE)* the upper 
surface trailing edge ordinate, since the last surface update. Should go to zero if 
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converging. In the viscous analysis case, it is the current value of YU(ITE); and 
it should approach a constant value if the solution is converging. Only changes 
after the first fifty cycles on each grid. 

SSP AT X — The current x ordinate value for the upper surface separation 
point ( x ■ -0.5 is the leading edge and x ■ 0.5 is the trailing edge). 

9. Boundary Layer Information 

Every LP cycles results of the current boundary layer solution are printed 
first for the lower surface of the airfoil and then for the upper surface. In each 
casei the laminar solution ( if HAM*!) is printed first followed by the turbulent 
solution. 

9a. Laminar Boundary Layer Information (Printed every LP cycles) 

X — Horizontal ordinate 
MACH # — Local Mach number 

CF — Skin friction coefficient (the O.iEii initial value is arbitrary and 
should be ignored) 

D-STAR — S /Ci non-dimensional boundary layer displacement thickness 
D-THETA — 0/c, non-dimensional momentum thickness 
H — 6*/Q shape factor 

RE-THETA — Local Reynolds number based on 0 
RE-STAR — Local Reynolds number based on 6* 

TM — 0* (du/dx)/T/ , Pressure gradient parameter 

9b. Possible Laminar Flow Messages 

a. SEPARATION OCCURRED AT X « O.xxx — gives x location where laminar 
separation occurred. 
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b. SHORT BUBBLE FORMED? TRANSITION TO TURBULENT FLOW 
ASSUMED, X- O.xxx 

c. LONO BUBBLE? LAMINAR STALL MAY OCCUR, X - O.xxx 
BOUNDARY LAYER CALCULATION WILL BE CONTINUED AS 
TURBULENT BUT ACCURACY OF RESULTS IS QUESTIONABLE 

d. BOUNDARY LAYER CALCULATION COMPLETED? 

NEITHER SEPARATION NOR TRANSITION WAS DETECTED 

9c. Turbulent Boundary Layer Information (Printed every LP cycles) 

X — Horizontal ordinate (-0.5 is leading edge, 0.5 is trailing edge) 

MACH # — Local Mach number 

DLSTR — 5*/c, non-dimensional boundary layer displacement thickness 
DEL — 6/c , non-dimensional boundary layer thicknesses 

CUE U e> -transformed boundary layer edge velocity, U 8 ^/^ >u 
USTAR — Law of the wall parameter, SQRT(tv/f ) 

USTAR**2 — Skin friction parameter, 

THETA — 0/c, ncn-dimensional momentum thickness 
CDP — Profile drag coefficient using Squire-Young approach modified 
according to Reference (9) 

CDW — Wave drag coefficient, not computed in this version of code and thus 
always zero 

CDTOT — Toxal drag coefficient 

9d. Turublent Boundary Layer Messages 

a. SEPARATED CP IS O.xxx — In the constant separated pressure level 
option (IMASS=l,ITACT 3 l,KSEP=i) the pressure m the separated region is printed 
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After each boundary layer calculation. If solution is converging* it should 
appraoch a constant value. 

b. USTAR2 (or USTCK) LT ZERO — Indicates where in computation separation 
was first detected. See program listing for details. 

10. Final Boundary Layer Results 

YUORIG — Original airfoil upper surface ordinate 
DU — Smoothed upper surface displacement thickness 
SLU — Slope of upper displacement surface 
YLORIG — Original airfoil lower surface ordinates 
DL — Smoothed lower surface displacement thickness 
SLL — Slope of lower displacement surface 

11. Pressure Distribution on Airfoil 

CPU — Upper surface pressure coefficient 
CPL — Lower surface pressure coefficient 

12. Final Displacement Surface Information 

YU — Ordinate of upper displacement surface 
YL — Ordinate of lower displacement surface 
SLU — Slope of upper displacement surface 

SLL — Slope of lower displacement surface i 

13. Mach Chart 

The Mach numbers at the I,J coordinate points are multiplied by 100 and 
printed out in block form. The grid points "inside" the upper and lower 

! 
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displacement surfaces are indicated by 2 eros. Velocities (U»V) at the flowfield 
grid points may also be printed out using option lPRT2*i. 

14. Miscellaneous Information 

The normal force coefficient) CN» and the drag coefficient) WAVE CD# 
obtained by integration of the pressure distribution are printed out. The latter 
should theoretically be zero for subcritical unseparated cases> but it usually is 
non-zero due to mesh size and grid placement. Thus values of WAVE CO should 
only be used for comparison purposes. 

15. Printer Plot of Results 

U — Upper surface pressure coefficient 
L — Lower surface pressure coefficient 
T — Upper displacement surface 
B — Lower displacement surface 

CPSTAR — Pressure coefficient for local Mach number of one 
CLCIR — Lift coefficient from computed circulation 
CL — Lift coefficient from integi ation of pressure distribution 
CD — Wave drag plus profile drag coefficient (Accuracy depends upon value 
of CDWAVE) 

CMLE — Moment coefficient about the leading edge 

CDF — Profile drag coefficient using Squire-Young approach modified for 
separation and compressibility 

CMC4 — Moment Coef x icient about quarter chord point 
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16. Miscellaneous Messages 


PROD .LB. ZERO — This message indicates that the turbulent boundary 
calculation was unable to obtain an appropriate starting solution. As a result the 
flowfield calculations were continued with the displacement surface ordinates 
and slopes -frozen at their last updated values. This situation usually occurs 
when transition takes place at the shock wave and the local Mach number 
immediately upstream of the shock wave is greater than 1.35. It probably 
indicates that at a minimum there is local separation m the vicinity of the shock 
wave. Since this phenomena is not modeled m the present code the results 
obtained in such cases should be used carefully. Sometimes this problem can be 
"avoided" by computing an all turbulent case (ILAM s 0) with transition at or 
immediately upstream of the shock wave. It should be noted that when this 
message appears, the subsequent computed values of CDF are in error. 

AA .LT. ZERO — This message occurs when the local cpeed of sound is 
computed to be negative. It indicates that either the solution has become 
unstable or else the rotated scheme tried to use a point outside of the solution 
domain due to a supersonic point on the JMAX-l row. Usual solution is to 
increase the supersonic dampingi EPSS, and/or increase the stretchi n and /or 
size of the Y grid. When encountered, the computation is termin.ced and the 
current perturbation flowfield solution is printed for diagnostic purposes. 

17. Note that additional output can be obtained for either diagnostic or analysis 
purposes by removing the C's from various commented print statements in the 
program. 
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TYPICAL RESULTS | 
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I 

In th® development of the present method, results have been obtained for ( 

NACA 4412 and N/.CA 0012 airfoils for freestream Mach numbers up to 0.6, angles 
of attack up to 16.5 degrees, and Reynolds numbers between three and nine 
million. These conditions were selected not because of limitations in the method 
but due to the availability of excellent experimental pressure distribution data in 
those regions. 19-21 Thus, the results presented here are only meant to be 
representative. j 

t 

Figures 4 and 5 compare results obtained with the present method with the 
low speed experimental data of Pinkerton^ 1 for a NACA 4412 airfoil at 6.3 million 
Reynolds number. In Figure 4 the experimental data has been plotted using the 
angle of attack correction suggested in Reference 21. As can be seen, the theory i 

predicts slightly larger lift coefficients than the data at the lower angles of i 

i 

attack. Whether this difference is due to an underestimation of the angle of 
attack correction, as suggested in Reference 22, or a problem in the theoretical 
model is unknown. In any event, the theory and the data are in excellent 
agreement between ten and fifteen degrees; and the present model reasonably 
predicts the location of maximum lift at an angle of attack of about 16 to 17 
degrees. The theoretical model does, however, overpredict slightly (1.9 vs. 1.7) 
the maximum lift coefficient predicted by this forty-seven year old data. 

Figure 5 compares pressure distribution results obtained with the present 
method with the experimental data of Pinkerton at an angle of attack slightly 
below that corresponding to maximum lift. In this case, the corrected angle of 
attack was used in the computations, and the upper and lower surface boundary 

37 


4i 


layers were assumed to be initially laminar followed by natural transtion to 
turbulent flow. For this high lift case* the theory predictes that the lower 
surface remaines entirely laminar and that the upper surface transitions at one 
percent chord followed by separation at 74.9 percent chord. As can be seen on the 
figure, the predicted pressure distribution is in excellent agreement with the 
data; and the pressure coefficient in the separated zone is slightly negative and 
constant. Experience indicates that for low speed cases better results are 
usually obtained using the constant pressure option <KSEP=0) for the separated 
zone. For this case, the theoretical lift coefficient was 1.69 while the 
experimental value was 1.68. The predicted profile drag coefficient was 0.0200, 
which is in reasor. ble agreement with availatle measurer., nts.23 

As indicated previously, it is important for this type of method to include the 
effects of a laminar boundary layer. Figure 6 shows lowar surface 
laminar-turbulent transition point locations predicted by the present method for 
a NACA 0012 airfoil at Mach 0.3 and six million Reynolds number. As can be seen, 
for these flight conditions the lower surface boundary layer is predominantly 
laminar at low angles of attack; and at angles of attack above ten degrees it is 
essentially all laminar. Obviously, a method which only includes a turbulent 
boundary layer calculation and/or assumes transition near the leading edge might 
yield incorrect results for these flight conditions. 

Figure 6 also shows for the same conditions the predicted upper surface 
separation points determined by the TAMSEP method. Notice that no upper 
surface separation is predicted until about 12 degrees angle of attack. After 
that, as the angle of attack increases, the beqinmng of separation moves forward 
on the upper surface until more than half of the airfoil experiences separated 
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flow at about 16 degrees. For this case* maximum lift occurs around fourteen 
degrees. 

Figure 7 compares the predicted lift coefficient variation with angle of attack 
for a NACA 0012 airfoil at the same nominal conditions (i.e. Hach 0.3 and six 
million Reynolds number) with experimental data obtained in the Low Turbulence 
Pressure Tunnel at NASA Langley.20 These data were obtained on clean airfoils 
without the use of trip strips* and thus the theoretical results were obtained 
using the laminar natural-transition turbulent model <ILAM a l). At the lower 
angles of attack the agreement is quite good. However* at tne onset of trailing 
edge separation* around twelve degrees* the theoretical lift curve exhibits a 
"kink" accompanied by a slight decrease in lift coefficient. This "kink" is often 
observed m the theoretical lift curves when separation is first predicted and is 
due to the fact that separation is usually first detected on the coarse grid. On 
that grid* the first point forward of the trailing edge is about 93% chord; r.nd* 
thus* when separation is first predicted on the coarse grid the separation point 
moves forward from the trailing edge to at least the 93 percent point* witn the 
result that the amount of separation is overpredicted and the lift at that 
condition is underpredicted. Since the model only permits the separation point to 
move forward* this effect is maintained thoughout all the grids at that angle of 
attack. However, as the angle of attack is increased this effect disappears. 
Thus, at low and medium freestream Mach numbers* the lift coefficients predicted 
at angles of attack just above the onset of trailing edge separation are usually 
slightly low. 

However, as can be seen on Figure 7* the lift coefficients predicted at the 
higher angles of attack are* at least for this case, in reasonable agreement v -ith 



th» experimental data. Notice that for this case tht experimental data indicates 

an apparent maximum lift around fourteen degrees followed by a dacraase and 
then an increata. Tha thaoretical li-ft basad upon the calculatad circulation 
pradicts a maximum lift coefficient of about 1.40 at 14 to Id degrees. On tha 
other hand* tha theoretical lift based upon integration of the pressure 
distributions indicates a maximum lift at fourteen degrees* which nominally 
agrees with experimental data. This slight divergence between the values 
predicted by circulation and those by pressure integration is* based upon 
experience* an excellent indication of the maximum lift location at low and 
medium speeds. The reason for this statement will be evident when the pressure 
distributions for these cases are discussed. 

The predicted variation of drag with lift fer 'his case is compared to 
experimental data at the higher lift coefficients on Figure 6. Again the 2 ig-zag 
in the theoretical curve corresponds to the similar phenomena on the lift versus 
angle of attacK plot and is du n to tne over;- .action of the size of the initial 
separated zone. Nevertheless* the agreement between the theoretical predictions 
and experimental values* particularly near maximum lift* is good and should be 
acceptable for applied engineering calculations. 

Figures 9<a-c) compare pressure distributions obtained with the present 
method with data obtained m the Low Turbulence Pressure Tunnel at If ASA 
Langley20 at three different ang);*? of attacK. The first corresponds to an 
unseparated flow situation* the second is near maximum lift* and the third is for 
an angle of attacK above the maximum lift condition. Since at medium freestream 
Mach numbers and hiqher evidence 9 ' 19 " 20 indicates that the presiure variation in 
a separated zone is not constant, these and subsequent cases were all run using 
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thi variable pressure option (KSEP«i). As can be seen on Figure 9(a) the 
theoretical pressure distribution -for the unsepsrsted case is in excellent 
agreement with experimental data. For this case* the theoretical lift and drag 
coefficients were 1.27 and 0.0116 while the corresponding clean airfoil 
experimental values were 1.23 and 0.0123. 

For the ca»e near maximum lift* Figure 9(b), it should be noted that some 
supersonic flow exists over the upper surface of the airfoil in a very small region 
near the leading edge since the critical C p tor this case is -6.09. In addition, the 
theoretical method predicts upper surface separation at 74.9 percent chord and 
boundary layer instability on the lower surface at about 60 percent chord. 
However, due to the favorable pressure gradient, the lower surface boundary 
layer never transitions. For this case the theoretical lift coefficient of 1.39 
coincides with the experimentally measured value, and the two pressure 
distributions exhibit reasonable agreement. 

At an angle of attack greater than that corresponding to maximum lift, the 
flow about an airfoil is typically characterized by a large region of unsteady 
separated flow; and a steady state solution method such as the present one is not 
really applicable. Thus, the apparent lack of agreement between the present 
steady theory and the experimental measurements shown on Figure 9tc» is not 
surprising. Nevertheless, the general pattern of the pressure distribution 
including the existence of a large separation zone is predicted; and the predicted 
lift coefficient of i .44 based upon circulation is in surprising agreement with the 
experimentally measured value of 1.437. However, careful examination of the 
solution indicates that it is not completely converged, that the theoretical lift 
may be oscillating slightly, and that the lift based upon pressure integration 
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computed 4t th» end of a run using default parameters is only 1.27. Interestingly* 
results using the thin layer Havier Stokes »quationsi7il8 ^ or 4 similar case 
(NACA 0012* Mach No. ® 0.3* Reynolds » 1 million* angle of attack ■ 18 degrees) 
indicate a lift coefficient varying with time from 0.65 to 1.6 with a Strouhal 
number of 0.1; and both the theoretical and experimental pressure profiles shown 
on Figure 9(c) are representative of those computed at various times in the cycle 
with the th n layer Navier Stokes model. Thus* the present theoretical result is 
representative of the type of pressure distribution and lift which might exist for 
this condition. 

Another interesting feature associated with the results shown on Figure 9(c) 
is that tht drag coefficient predicted using the methrd of Reference 6 (i.e. CDF) 
was 0.0601 while that predicted using the method of Reference 9 ( i.e. CDP) was 
0.1190. Kormally* these two values are m good agreement with each other. 
Apparently* in the present method when the maximum lift condition is exceeded 
the solution becomes oscillatory and not completely converged* the lift computed 
by pressure integration diverges from and is lower than that from circulation* and 
the CDF and CDP values differ significantly. It is believed that these three 
items can be used to determine for medium Mach numbers the angle of attack 
corresponding to maximum lift. 

Figure 10 shows for the same NACA 0012* Mach 0.3, Reynolds number 6 
million case, the displacement surfaces predicted by the present method for the 
upper surface region between 70 percent chord and the trailing edge at various 
angles of attack. Below 11.13 degrees, where the flow is unseparated* the 
displacement thicknesses in the trailing edge zone are relatively small. However, 
with the onset of separation at 12.09 degrees the thicknesses begin to increase 
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rapidly* and the displacement surfaces take on shapes characteristic of the flow 
over a stalled airfoil. It should be noticed that at an angle of attack of id. 12 
degrees* the displacement surface starts to curve back towards the freestream 
angle prior to the trailing edge. It is believed that the displacement surfaces 
shown on this figure have the correct behavior and are an adequate engineering 
representation of the real flow. 

f 

♦ 

l 

Obviously* it would be desirable for the present method to accurately model 
transonic flows with and without significant trailing edge separation. 
Consequently* predictions obtained with TAM SEP have been compared with data 
| obtained in the NASA 8 Ft. Transonic Pressure Tunnel by Harris. 19 *j>his j&ta 

probably represents the best high lift transonic experimental airfoil data 
available today* and it has been used by many investigators. However* in most 
cases previous studies have compared results using the tunnel geometric angle of 
attack values and have ignored the known corrections associated with this data. 
Since such comparisons could lead to erroneous conclusions* the present results 
were obtained by matching the tunnel Mach and Reynolds numbers and using the 
corrected angles of attack suggested by Harris. 

Figure 11 shows a transonic separated flow result compared with data 
obtained by Harris. 19 while the pressure distribution shown was obtained using 
the laminar-turbulent boundary layer option* indistinguishable results were 
‘ obtained assuming transition at six percent chord. In the actual experiment* 

| boundary layer trip strips were located at five percent chord. For this case* the 

* 

I 

. present method predicts upper surface separation at 87 percent chord and lower 

* surface transition very near the trailing edge. While the experimental shock 

! location is slightly forward of the theoretical value and while the experimental 
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pressures in the trailing edge region ere slightly lower then those predicted by 
the theory, the overell egreement is probebly eccepteble tor engineering studies. 
For this case, the measured normel force coefficient we* 0.994 while the 
predicted lift coefficient wes 0.981. 

A lift versus engle of attack curve typical of those predicted by the present 
method is compered with experimentel dete on Figure 12. For this freestreem 
Mech number of 0.5i significant transonic flow accompanied by a strong shock 
wave is present on the upper surface at angles of attack greater than six 
degrees. As can be seen, the theoretical prediction* which was obtained assuming 
transition at six percent chord* agrees well with the experimental data up to 
about 7.35 degrees. Above that angle of attack, the present method predicts a 
maximum lift coefficient of 1.09 at 8.24 degrees with no trailing edge separation. 
At 9.21 degrees the theory predicts a decrease in lift coefficient to C.982 with 

upper surface separation at 87 percent chord. The experimental data* however* 
indicates that maximum lift is i.02 at 9.21 degrees and that trailing edge 
separation probably started at about 7.S degrees. Examination of the theoretical 
results at 7.35 degrees reveals that the local Hach number immediately upstream 
of the upper surface shock wave is 1.42. Such a strong shock wave* whose 
strength increases with angle of attack, should induce significant shock boundary 
layer interaction which, unfortunately, is not modeled in the present theory and 
code. Therefore, in this case the differences between the theoretical rc-.ults and 
the experimental data near maximum lift are probably due to shock boundary layer 
interaction and its subsequent effect on boundary layer growth and trailing edge 
separation. Nevertheless, the present model and code does give a reasonable 
indication of the location and magnitude of the maximum lift coefficient. 
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A similar lift curve for a free stream Mach number of 0.55 is shown on Figure 
13. As before> the theoretical results were obtained using an all turbulent 
boundary layer; and in this case no upper surface trailing edge separation was 
detected until an angle of attack of 8.266 degrees. The maximum lift coefficient 
was computed to be 1.02 at 7.29 degrees as compared to the experimental values 
of 0.983 and 8.266 degrees. For this case significant transonic flow existed at 
all angles of attack above 4.5 degrees and by 7.29 degrees the Mach number at the 
upper surface shock wave had increased to 1.50. At 9.33 degrees a complete 
solution could not be obtained with the present method. On the medium grid the 
upper surface flow separated at 87 percent chord* and on the fine grid the 
separation point moved forward to the shock wave at about 18 percent chord and 
the solution failed. Quite obviously significant shock boundary layer interaction 
exists at these high angles of attack and the decrease in lift or stall is probably 
due more to shock induced separation than to the onset of significant trailing 
edge separation. 

Nevertheless* the present method can be used to estimate reasonably 
accurately the occurence of this situation. As can be seen on Figure 13* the 
method predicts reasonably well the magnitude of the maximum lift coefficient 
and is conservative as to the corresponding angle of attack location. In addition, 
by noting the mechanism of code "failure", in this case sudden separation at the 
shock wave* a user can probably determine the type of stall phenomena present. 

Theoretical pressure distribution results are compared with experimental data 
for a freestream Mach number 0.6 case m Figure 14. This case at a corrected 
angle of attack of 5.59 degrees and three million Reynolds number is significant 
for a variety of reasons. First, it is an example of a transonic case with a strong 
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upper surface shock wave. Second) the flow is unseparated and the data should 
serve as a good test of the present method for a situation without separation. 

Finally) this case has also been solved by Anderson et all7 US jng both an Buler 
boundary layer method and a thin layer Navier Stokes method. 

1 

! 

It should be noticed that the present results) like those of Anderson> agree 
very well with the experimental pressure distribution with respect to shock 
location and pressure levels. Also) for this case the experimental lift coefficient 
was 0.78i. The present TAMSBP method predicted 0.80?) the Euler boundary layer 
method of Anderson yielded 0.804) while his thin layer Nav.er Stokes result was 
0.793. Obviously) the present method is capable of yield' ig excellent results that 
are in agreement with experimental data and other analytical methods. However) 
since it is a non-conservative full potential method) it should obtain such results 

i 

with accurate shock wave locations more easily and faster than other more 

i i 

complicated methods. 

Another lift versus angle of attack comparison is shown on Figure 15 for the ' 

t 

NACA 0012 at a freestream Mach number of 0.6 and a Reynolds number of nine 

1 

million. As can be seen the highest theoretical point plotted is at a lift 
coefficient value of 0.945 and an angle of attack of 6.392 degrees. At 7.348 
degrees the boundary layer solution became frozen (i.e. PROD.LE.2BRO) on the 
medium grid due to shock boundary layer interaction, and on the fine grid the flow 
separated at the shock wave and a converged solution was not obtained. At 6.371 
degrees the flow separated at the upper surface shock wave on the coarse grid, 
and again a converged solution was not obtained. For these cases, computed local 
Mach numbers in the vicinity of the shock wave were as high as 1.56. Thus, the 
theory indicates that above 6.392 degrees significant shock boundary layer 
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interaction probably accompanied by Reparation exists and that the maximum lift 
coefficient occurs at that angle of attacK. As can be seen on the figure* the 
magnitude of the predicted maximum lift coefficient is in good agreement with the 
experimental data; although the angle of attack location is again conservative. 

Based upon the results presented in this section* it is believed that the 
present method and code can be used at low and medium Mach numbers to 
accurately predict lift and pressure distributions at angles of attack up to that 
associated with maximum lift. At transonic speeds* tho method should give good 
results for unseparated flows and for flows having trailing edge separation 
without significant shock boundary layer interaction. Thus* at transonic 
conditions* the method is probably currently limited* for accurate results* to 
Reynolds numbers of three million and higher and to local upper surface Mach 
numbers less than 1.4 to 1.45. In addition* it should yield reasonable estimates 
for the maximum lift coefficient at transonic speeds while being conservative as 
to the corresponding value of angle of attack; and the method should indicate the 
onset of significant shack boundary layer interaction. 
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CONCLUSION 


A direct-invars* technique based upon a nonconservative full potential 
inviscid method* a Thwaitss laminar boundary layer technique* and the Barnwell 
turbulent momentum integral method has been developed. This method is suitable 
for predicting the subsonic and transonic flowfield about airfoils having trailing 
edge separated flow. Extensive comparisons with experimental data indicate that 
the method should be a useful tool for applied aerodynamic engineering analyses. 
In addition* it is believed that the range of applicability of the method could be 
extended significantly by the addition of a shock boundary layer interaction mode. 
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Figure 1 — Problem Formulation 
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Figure 2 — Iteration and Interaction Procedure Used in TAMSEP 
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Figure 3 — Subroutine Tree for TAMSEP 
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Figure 4 — Comparison of Predicted Lift Coefficient with Experimental Data 

NACA 4412 Case 
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Figure 7 — Comparison of Predicted Lift Coefficient with Experimental Data 
NACA 0012 Mach No. 3 0.3 Case 
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Figure 8 — Comparison of Theoretical Lift and Drag with Experimental Data 
NACA 0012 Mach No. » 0.3 Case 
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Figure 9(b) — Theoretical and Experimental Pressure Distribution Comparisons 
NACA 0012 at 14.3 Degrees Angle of Attack 
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Figure 14 -- Theoretical and Experimental Pressure Distribution Comparison 

NACA 0012 at Mach No. 0.6 
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